#ifndef SPLINTER_PPFORM_H
#define SPLINTER_PPFORM_H

#include "function.h"

#include "polynomial.h"

namespace SPLINTER
{

    /**
     * Class that implements the multivariate tensor product ppForm
     */
    class SPLINTER_API PPForm : public Function
    {
    public:
        /**
         * Builder class for construction by regression
         * Implemented in ppFormBuilder.*
         */
        class Builder;
        enum class Smoothing;
        enum class KnotSpacing;

        PPForm(unsigned int numVariables);

        /**
         * Construct ppForm from knot vectors, coefficients, and basis degrees
         */
        PPForm(std::vector<std::vector<double>> knotVectors, std::vector<unsigned int> basisDegrees);
        PPForm(std::vector<std::vector<double>> polynomials, std::vector<std::vector<double>> knotVectors, std::vector<unsigned int> basisDegrees);
        PPForm(std::vector<Polynomial> polynomials, std::vector<std::vector<double>> knotVectors, std::vector<unsigned int> basisDegrees);

        /**
         * Construct ppForm from file
         */
        PPForm(const char *fileName);
        PPForm(const std::string &fileName);

        virtual PPForm *clone() const { return new PPForm(*this); }

        /**
         * Evaluation of ppForm
         */

        // Avoid name hiding
        using Function::eval;
        using Function::evalHessian;
        using Function::evalJacobian;

        // Evaluation of ppForm
        double eval(DenseVector x) const override;
        DenseMatrix evalJacobian(DenseVector x) const override;
        DenseMatrix evalHessian(DenseVector x) const override;

        // Evaluation of ppForm basis functions
        SparseVector evalBasis(DenseVector x) const;
        // SparseMatrix evalBasisJacobian(DenseVector x) const;

        /**
         * Getters ;
         */
        std::vector<Polynomial> getPolynomials()
        {
            return polynomials;
        }

        unsigned int getNumPolynomials() const
        {
            return polynomials.size();
        }

        unsigned int getNumControlPoints() const
        {
            return polynomials.size();
        }

        // std::vector<unsigned int> getNumBasisFunctionsPerVariable() const;

        DenseMatrix getControlPoints() const;
        std::vector<std::vector<double>> getKnotVectors() const;
        std::vector<unsigned int> getBasisDegrees() const { return basisDegrees; }
        // std::vector<double> getDomainUpperBound() const;
        // std::vector<double> getDomainLowerBound() const;

        /**
         * Setters
         */
        void setPolynomials(const std::vector<Polynomial> &polynomials);
        void setControlPoints(const DenseMatrix &controlPoints);
        void checkControlPoints() const;

        // Linear transformation of control points (ppForm has affine invariance)
        void updateControlPoints(const DenseMatrix &A);

        // Reduce support of ppForm
        // void reduceSupport(std::vector<double> lb, std::vector<double> ub, bool doRegularizeKnotVectors = true);

        // Perform global knot refinement
        // void globalKnotRefinement(); // All knots in one shabang

        // Perform a local knot refinement at x
        // void localKnotRefinement(DenseVector x);

        // Decompose ppForm to Bezier form
        // void decomposeToBezierForm();

        // Insert a knot until desired knot multiplicity is obtained
        // void insertKnots(double tau, unsigned int dim, unsigned int multiplicity = 1);

        void save(const std::string &fileName) const override;

        std::string getDescription() const override;

    protected:
        PPForm();

        /*
         * The control point matrix is P = (knotaverages, coefficients) in R^(m x n),
         * where m = numBasisFunctions and n = numVariables + 1. Each row in P is a control point.
         */
        std::vector<Polynomial> polynomials;
        DenseMatrix knotaverages;
        std::vector<std::vector<double>> knotVectors;
        std::vector<unsigned int> basisDegrees;

        // Control point computations
        // DenseMatrix computeKnotAverages() const;

    private:
        // Domain reduction
        // void regularizeKnotVectors(std::vector<double> &lb, std::vector<double> &ub);
        // bool removeUnsupportedBasisFunctions(std::vector<double> &lb, std::vector<double> &ub);

        // Helper functions
        // bool pointInDomain(DenseVector x) const;

        void load(const std::string &fileName) override;

        friend class Serializer;
        friend bool operator==(const PPForm &lhs, const PPForm &rhs);
    };

} // namespace SPLINTER

#endif // SPLINTER_PPFORM_H